################################################################################
############################07. Figures 2 and 3#################################
################################################################################
library(MASS)
library(gridExtra)

#Make dataset for the rug plots

pols_rug <- readRDS('pols_bjps') %>%
  dplyr::select(scale_econsdw,scale_dispgini, 
                scale_ggini_inc1,
                scale_socioeconomicsal,
                scale_galtansdw ,
                scale_enep ,
                scale_unemp) %>%
  na.omit
#####Figure 2: Sorting on X-axis, salience held constant in each panel#####

#simulate 1000 draws from a multivariate random normal distribution based on the model#
test <- as.data.frame(mvrnorm(1000, mu =coef(full1), Sigma =vcovHC(full1)))

#vector of sorting from observed min to observed max#
nums <- seq(from = min(pols$scale_ggini_inc1, na.rm = T), to = max(pols$scale_ggini_inc1, na.rm = T), length.out = 1000)

#create data frame for storing output#
marg <- data.frame(x = nums)

#Simulate marginal effect based on draws from the model#
for(i in 1:1000){
  for(j in 1:1000){
    marg[i,j+1] <- test[j,1] + test[j,7]*marg[i,1] + test[j,8]*-1 + test[j,10]*marg[i,1]*-1
  }
}

#Vectors to be filled#
sort <- NA
pe <- NA
lb <- NA
ub <- NA


#Make CI based on empirical quantiles

for(i in 1:1000){
  sort[i] <- marg$x[i]
  pe[i] <- unlist(quantile(marg[i, 2:1001], probs = .5))
  lb[i] <- unlist(quantile(marg[i, 2:1001], probs = .025))
  ub[i] <- unlist(quantile(marg[i, 2:1001], probs = .975))
}

#Combine into a singe dataset

marg.plot.dat <- data.frame(sort = nums,
                            pe = unlist(pe),
                            lb = unlist(lb),
                            ub = unlist(ub))

#create plot. 

lwr.plot.v <- ggplot(marg.plot.dat, aes(x = sort, y = pe)) +
  geom_ribbon(aes(ymin = lb, ymax = ub), fill = 'grey69', alpha = .5) +
  geom_line() +
  ylim(c(-1.5, 3.5))+
  xlim(c(-2,2))+
  geom_hline(aes(yintercept = 0), linetype = 'longdash', color = 'red') +
  labs(title = "Economic Salience = -1", x = "", y = "") +
  theme_bw() +
  geom_rug(data = pols_rug, aes(x = scale_ggini_inc1), inherit.aes = F)

#These steps are repeated for each of the panels in each figure.#

#####Middle Panel####

for(i in 1:1000){
  for(j in 1:1000){
    marg[i,j+1] <- test[j,1] + test[j,7]*marg[i,1] + test[j,8]*0 + test[j,10]*marg[i,1]*0
  }
}

sort <- NA
pe <- NA
lb <- NA
ub <- NA



for(i in 1:1000){
  sort[i] <- marg$x[i]
  pe[i] <- unlist(quantile(marg[i, 2:1001], probs = .5))
  lb[i] <- unlist(quantile(marg[i, 2:1001], probs = .025))
  ub[i] <- unlist(quantile(marg[i, 2:1001], probs = .975))
}

marg.plot.dat <- data.frame(sort = nums,
                            pe = unlist(pe),
                            lb = unlist(lb),
                            ub = unlist(ub))

middle.plot.v <- ggplot(marg.plot.dat, aes(x = sort, y = pe)) +
  geom_ribbon(aes(ymin = lb, ymax = ub), fill = 'grey69', alpha = .5) +
  geom_line() +
  ylim(c(-1.5, 3.5))+
  xlim(c(-2,2))+
  geom_hline(aes(yintercept = 0), linetype = 'longdash', color = 'red') +
  labs(title = "Economic Salience = 0", x = "", y = "") +
  theme_bw()+
  geom_rug(data = pols_rug, aes(x = scale_ggini_inc1), inherit.aes = F)

####Right Plot#####

for(i in 1:1000){
  for(j in 1:1000){
    marg[i,j+1] <- test[j,1] + test[j,7]*marg[i,1] + test[j,8]*1 + test[j,10]*marg[i,1]*1
  }
}

sort <- NA
pe <- NA
lb <- NA
ub <- NA



for(i in 1:1000){
  sort[i] <- marg$x[i]
  pe[i] <- unlist(quantile(marg[i, 2:1001], probs = .5))
  lb[i] <- unlist(quantile(marg[i, 2:1001], probs = .025))
  ub[i] <- unlist(quantile(marg[i, 2:1001], probs = .975))
}

marg.plot.dat <- data.frame(sort = nums,
                            pe = unlist(pe),
                            lb = unlist(lb),
                            ub = unlist(ub))

righ.plot.v <- ggplot(marg.plot.dat, aes(x = sort, y = pe)) +
  geom_ribbon(aes(ymin = lb, ymax = ub), fill = 'grey69', alpha = .5) +
  geom_line() +
  ylim(c(-1.5, 3.5))+
  xlim(c(-2,2))+
  geom_hline(aes(yintercept = 0), linetype = 'longdash', color = 'red') +
  labs(title = "Economic Salience = 1", x = "", y = "") +
  theme_bw()+
  geom_rug(data = pols_rug, aes(x = scale_ggini_inc1), inherit.aes = F)

Figure_3_fe <- grid.arrange(lwr.plot.v, middle.plot.v, righ.plot.v, nrow = 1, top = "", left = "Marginal Effect of Disposable Income Inequality", bottom = "Sorting by Income")

#ggsave("figure_3_fe.pdf", plot = Figure_2_fe, width = 11, height = 5, units = 'in')

#####Figure 3: Salience on x-axis, sorting held constant#####


#vector of sorting from observed min to observed max#
nums <- seq(from = min(pols$scale_socioeconomicsal, na.rm = T), to = max(pols$scale_socioeconomicsal, na.rm = T), length.out = 1000)

#create data frame for storing output#
marg <- data.frame(x = nums)

#Simulate marginal effect based on draws from the model#
for(i in 1:1000){
  for(j in 1:1000){
    marg[i,j+1] <- test[j,1] + test[j,7]*-1 + test[j,8]*marg[i,1] + test[j,10]*marg[i,1]*-1
  }
}

#Vectors to be filled#
sal <- NA
pe <- NA
lb <- NA
ub <- NA


#Make CI based on empirical quantiles

for(i in 1:1000){
  sal[i] <- marg$x[i]
  pe[i] <- unlist(quantile(marg[i, 2:1001], probs = .5))
  lb[i] <- unlist(quantile(marg[i, 2:1001], probs = .025))
  ub[i] <- unlist(quantile(marg[i, 2:1001], probs = .975))
}

#Combine into a singe dataset

marg.plot.dat <- data.frame(sal = nums,
                            pe = unlist(pe),
                            lb = unlist(lb),
                            ub = unlist(ub))

#create plot. 

lwr.plot.v <- ggplot(marg.plot.dat, aes(x = sal, y = pe)) +
  geom_ribbon(aes(ymin = lb, ymax = ub), fill = 'grey69', alpha = .5) +
  geom_line() +
  ylim(c(-1.5, 3.5))+
  xlim(c(-2,2))+
  geom_hline(aes(yintercept = 0), linetype = 'longdash', color = 'red') +
  labs(title = "Sorting = -1", x = "", y = "") +
  theme_bw()+
  geom_rug(data = pols_rug, aes(x = scale_socioeconomicsal), inherit.aes = F)

#These steps are repeated for each of the panels in each figure.#

#####Middle Panel####

for(i in 1:1000){
  for(j in 1:1000){
    marg[i,j+1] <- test[j,1] + test[j,7]*0 + test[j,8]*marg[i,1] + test[j,10]*marg[i,1]*0
  }
}

sal <- NA
pe <- NA
lb <- NA
ub <- NA



for(i in 1:1000){
  sort[i] <- marg$x[i]
  pe[i] <- unlist(quantile(marg[i, 2:1001], probs = .5))
  lb[i] <- unlist(quantile(marg[i, 2:1001], probs = .025))
  ub[i] <- unlist(quantile(marg[i, 2:1001], probs = .975))
}

marg.plot.dat <- data.frame(sal = nums,
                            pe = unlist(pe),
                            lb = unlist(lb),
                            ub = unlist(ub))

middle.plot.v <- ggplot(marg.plot.dat, aes(x = sal, y = pe)) +
  geom_ribbon(aes(ymin = lb, ymax = ub), fill = 'grey69', alpha = .5) +
  geom_line() +
  ylim(c(-1.5, 3.5))+
  xlim(c(-2,2))+
  geom_hline(aes(yintercept = 0), linetype = 'longdash', color = 'red') +
  labs(title = "Sorting = 0", x = "", y = "") +
  theme_bw()+
  geom_rug(data = pols_rug, aes(x = scale_socioeconomicsal), inherit.aes = F)

####Right Plot#####

for(i in 1:1000){
  for(j in 1:1000){
    marg[i,j+1] <- test[j,1] + test[j,7]*1 + test[j,8]*marg[i,1] + test[j,10]*marg[i,1]*1
  }
}

sal <- NA
pe <- NA
lb <- NA
ub <- NA



for(i in 1:1000){
  sal[i] <- marg$x[i]
  pe[i] <- unlist(quantile(marg[i, 2:1001], probs = .5))
  lb[i] <- unlist(quantile(marg[i, 2:1001], probs = .025))
  ub[i] <- unlist(quantile(marg[i, 2:1001], probs = .975))
}

marg.plot.dat <- data.frame(al = nums,
                            pe = unlist(pe),
                            lb = unlist(lb),
                            ub = unlist(ub))

righ.plot.v <- ggplot(marg.plot.dat, aes(x = sal, y = pe)) +
  geom_ribbon(aes(ymin = lb, ymax = ub), fill = 'grey69', alpha = .5) +
  geom_line() +
  ylim(c(-1.5, 3.5))+
  xlim(c(-2,2))+
  geom_hline(aes(yintercept = 0), linetype = 'longdash', color = 'red') +
  labs(title = "Sorting = 1", x = "", y = "") +
  theme_bw()+
  geom_rug(data = pols_rug, aes(x = scale_socioeconomicsal), inherit.aes = F)

Figure_4_fe <- grid.arrange(lwr.plot.v, middle.plot.v, righ.plot.v, nrow = 1, top = "", left = "Marginal Effect of Disposable Income Inequality", bottom = "Salience of Economic Issues")

#ggsave("figure_4_fe.pdf", plot = Figure_3_fe, width = 11, height = 5, units = 'in')